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Abstract 

We present the details of an algorithm for the global evolution of asymptot- 
ically flat, axisymmetric spacetimes, based upon a characteristic initial value 
formulation using null cones as evolution hypersurfaces. We identify a new 
static solution of the vacuum field equations which provides an important 
test bed for characteristic evolution codes. We also show how linearized solu- 
tions of the Bondi equations can be generated by solutions of the scalar wave 
equation, thus providing a complete set of test beds in the weak field regime. 
These tools are used to establish that the algorithm is second order accurate 
and stable, subject to a Courant-Friedrichs-Lewy condition. In addition, the 
numerical versions of the Bondi mass and news function, calculated at scri 
on a compactified grid, are shown to satisfy the Bondi mass loss equation 
to second order accuracy. This verifies that numerical evolution preserves 
the Bianchi identities. Results of numerical evolution confirm the theorem 
of Christodoulou and Klainerman that in vacuum, weak initial data evolve 
to a flat spacetime. For the class of asymptotically flat, axisymmetric vac- 
uum spacetimes, for which no nonsingular analytic solutions are known, the 
algorithm provides highly accurate solutions throughout the regime in which 
neither caustics nor horizons form. 
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I. INTRODUCTION 



The physical basis of a new algorithm for the evolution of spacetimes has been pro- 
posed. |IJ This algorithm is based upon the characteristic initial value problem for general 
relativity, using light cones as the evolution hypersurfaces, rather than the spacelike folia- 
tion used in traditional approaches based upon the Cauchy problem. The intimate use of 
characteristics has particular advantages for the description of gravitational radiation and 
black hole formation. [0-H|. The first attempt to carry out numerical evolution based upon 
this algorithm was only successful in a region outside some sufficiently large worldtube. Near 
the vertex of the null cone, instabilities arose which destroyed the accuracy of the code. The 
underlying cause of this instability was too complicated to analyze in the context of general 
relativity, especially since the numerical analysis of the characteristic initial value problem 
had not yet been carried out even for the simplest linear axisymmetric systems. 

This warranted an investigation of the basic computational properties of the evolution 
of the flat space scalar wave equation using a null cone initial value formulation It was 
found that near the vertex of the cone the Courant-Friedrichs-Lewy (CFL) condition places 
a stricter limit on the time step than for the case of Cauchy evolution. This insight made 
possible the development of an extremely efficient marching algorithm for evolving the data 
on the initial cone by stepping it out from the vertex of the next cone to null infinity (scri) 
along each angular ray direction. This marching algorithm is based upon a simple identity 
satisfied by the values of the scalar field at the corners of a parallelogram formed by four 
null rays. The result was a stable, calibrated, second order accurate global algorithm on 
a compactified grid. Furthermore, scri played the role of a perfect absorbing boundary so 
that no radiation was reflected back into the system. This algorithm offers a powerful new 
approach to generic wave type systems. The basic idea is applicable to any of the hyperbolic 
systems occurring in physics. 

In this paper, we apply the algorithm to the evolution of asymptotically flat, (twist-free) 
axisymmetric spacetimes. In a Bondi null coordinate system the metric takes the form |J 



The vacuum field equations then decompose into the three hypersurface equations 



ds 2 = (-e 2f3 - UW^du 2 + 2e 2f3 dudr + 2Ur 2 e 2 ^dud9 



r 



-r 2 {e 2 ^dd 2 + e~^ sin 2 9d(j) 2 ). 



(1) 



P,r = \r ( 7 , r ) 2 



(2) 




(3) 




(4) 
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and one evolution equation 



4r(r 7 ) iUr = [2r 7>r V - r 2 (2 7 , e 17 + sin 0(__) „)] - 2r 

sin sm # 

+ l r V(T-«(cr, r ) a + 2e 2 ^[(A,) 2 + sin0(^U (5) 

The initial data consists of 7, which is unconstrained except by smoothness conditions. 
Because 7 represents a spin-2 field, it must be 0(sin 2 #) near the axis and consist of I > 2 
spin-2 multipoles. 

Here we are interested in the global application of this system when the null hypersur- 
faces are null cones, although the approach also goes through without major change if the 
null hypersurfaces emanate from a finite world tube. We require that the null cones have 
nonsingular vertices which trace out a geodesic worldline r = 0. For the quadrupole terms, 
this implies the boundary conditions 7 = 0(r 2 ), (3 = 0(r 4 ), U = 0(r) and V = r + 0(r 3 ). 
For higher multipoles, the smoothness conditions can be worked out by referring back to 
local Minkowski coordinates 0. As a consequence, 0(r n ) terms in 7 can contain multipoles 
with 2 < I < n. Any satisfactory computational algorithm must meet the challenge of 
preserving these smoothness conditions. 

In Sec. we analyze the linearized version of these equations and show how their 
solutions may be obtained locally from solutions of the scalar wave equation. This provides 
an important means of constructing linearized solutions in a null cone gauge for the purpose 
of code calibration. It also reveals essential changes in the grid structure necessary in 
adapting the null parallelogram algorithm for the wave equation to linearized gravity. 



This also solves the major computational problems for the nonlinear case. In Sec. Ill 



we show how the linearized algorithm can be extended to this case. In Sec. [IV], we discuss 
the major finite difference techniques necessary for second order accuracy. In Sec. [V], we 
present a study of the stability and accuracy of an evolution code based upon this algorithm. 
New exact and linearized solutions are introduced to establish second order convergence. In 
addition, a global check on accuracy is carried out using Bondi's formula relating mass loss 
to the time integral of the square of the news function. 

II. THE LINEARIZED BONDI EQUATIONS 

In the linearized limit (3 = and V — r. The equations ©-(HD reduce to one hypersurface 
equation and one evolution equation for the surviving field variables 7 and U, 

{r%rh = ( 6) 
sin 

4r(r 7 ) jOT = [2r 2 7jr - r 2 sin#(^-\ e ] ir . (7) 

sin u 

Physical considerations suggest that these equations be related to the wave equation. If 
this relationship were sufficiently simple, then the scalar wave algorithm could be used as a 
guide in formulating an algorithm for evolving 7. A scheme for generating solutions to the 
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linearized Bondi equations in terms of solutions to the wave equation has been presented 
previously ||. However, in that scheme, the relationship of the scalar wave to 7 is non-local 
in the angular directions and is not useful for this purpose. 

We now formulate an alternative scheme in terms of spin-weight quantities a and Z, 
related to 7 (spin-weight 2) and U (spin-weight 1) by [|J 

7 = 5 2 a = sm9d e (—^— d e a) (8) 
sin 9 



U = 3Z = d e Z. (9) 
Then the linearized equations are equivalent to 

(r 4 Z r ), r = -2r 2 (2 - L 2 )a r , (10) 

and 

E := 2(ra), ur - r _1 (r 2 a r - r 2 Z/2) >r = 0, (11) 

where L 2 = —(1/ sin 8)dg(sin 9dg) is the #-part of the angular momentum operator. Now let 
$ be a solution of the flat space scalar wave equation, 

rD $ = 2(r$) iUr - (r$) )JT + r _1 L 2 $ = 0, (12) 



and set 



and 



Then 



and 



r 2 cv = (r 2 $), r (13) 



r 2 Z iT . = 2(L 2 - 2)$. (14) 



E = rO $ + 2($ + ar),« - 2r~ 2 (r 2 $) ir + Z, (15) 



S = r- 2 (r 3 D$) jr (16) 



Equation (|T0|) is satisfied as a result of (]T5|) and (13) and the wave equation (|T2|) implies 
_E r = 0. If $ is smooth and 0{r 2 ) at the origin, this implies E = 0, so that the linearized 
equations are satisfied. The condition that $ = 0(r 2 ) eliminates fields with only monopole 
and dipole dependence so that it does not restrict the generality of the spin-weight 2 function 
7 obtained through (|). 

Thus for any linearized axisymmetric gravitational wave in the null cone gauge, 7 may 
be related to a scalar wave by and (|i~3"|). The CFL condition for convergence of a 
finite difference algorithm requires that the numerical domain of dependence be larger than 
the physical domain of dependence. Because the relationship between $ and 7 is local with 
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respect to the null rays on the cone, their domains of dependence coincide. This suggests that 
a stable and convergent evolution algorithm for the gravitational field may be modeled upon 
the scalar wave algorithm. This turns out to be the case although some subtle distinctions 
arise, as described below. 

An evolution algorithm for scalar waves has been formulated in terms of an integral 
identity for the values of the field at the corners of a null parallelogram lying on the (u, r) 
plane 0. The wave equation with source, □$ = S, can be reexpressed in the form 

□ = + rS, (17) 

where ip = r<3> and □ ^ is the 2-dimensional wave operator intrinsic to the (u, r) plane. 
Integration over the null parallelogram as depicted in Fig. [I] then leads to the integral 
equation 

1 f L 2 ip 

ipQ = + i>S - + o / dudr[ — + rS], (18) 

where P, Q, R and S are the corners and A the area of the parallelogram. 

This identity gives rise to an explicit marching algorithm for evolution. Let the null 
parallelogram span null cones at adjacent grid values uq and uq + Am, as shown in Fig. [I], 
for some 9 and <fi. If ip has been determined on the entire uq cone and on the uq + Am cone 
radially outward from the origin to the point P, then fll8|) determines at the next point 
Q in terms of an integral over A. This procedure can then be repeated to determine ip at 
the next radially outward point T in Fig. [I]. After completing this radial march to scri, the 
field ip is then evaluated on the next null cone at uq + 2Am, beginning at the vertex where 
smoothness gives the start up condition that ip = 0. The resulting evolution algorithm is 
a 2-level scheme which reflects, in a natural way, the distinction between characteristic and 
Cauchy evolution, i.e. that the time derivative of the field is not part of the characteristic 
initial data. 

The CFL condition on the numerical domain of dependence is a necessary condition 
for convergence of a numerical algorithm. For the grid point at (u,r,9), there are three 
critical grid points, (u — An, r + Ar, 9) and (u — Au, r — Ar, 9 ± A9), which must lie inside 
its past physical light cone. These gives rise to the inequalities Au < 2Ar and Au < 
—Ar + (Ar 2 + r 2 A9 2 ) 1 ' 2 . At large r, the second inequality becomes Am < rA9 and the 
limitations on the time step are essentially the same as for a Cauchy evolution algorithm. 
However, near the vertex of the cone, the second inequality gives a stricter condition 

Am < KArA9 2 , (19) 

where K is a number of order 1 whose exact value depends upon the start up details at 
the vertex. For the scalar wave equation, these stability limits were confirmed by numerical 
experimentation and it was found that K w 4. 

The linearized gravitational evolution equation ([7|) can be recast into a form similar to 

(0), 

□ (2 V = H, (20) 
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where now ijj = r'y and Ti only contains hypersurface terms. This allows use of the null 
parallelogram algorithm to evolve 7 by the same marching scheme as in the scalar case. The 
additional feature here is that U must be simultaneously marched out the null cone using 
the hypersurface equation (|). For the scalar wave equation (|17D , the angular momentum 
barrier is represented by the term L 2 ip/r 2 , which is determined from ip by a double angular 
derivative. In the linearized gravitational evolution equation (0), the analogous term is 
[r 2 sin6(U/ sin 0)^}^, which is determined from U by a single angular derivative. In turn, 
the hypersurface equation (0) relates U to 7 by a single angular derivative. Physically, 
this has the same net effect of producing an angular momentum barrier depending upon 
the second angular derivative, as in the scalar case. However, the nontrivial mathematical 
distinction between the two cases leads to nontrivial difference in their natural grid structures 
for a numerical algorithm. In particular, the grid for U must be staggered between the grid 



points for 7. These and other details of the marching algorithm will be given in Sec. |TV], 
where we discuss the nonlinear case. 

The use of scalar waves to generate solutions of the linearized Bondi equations provides 
an important tool for testing evolution algorithms in the linear regime. Monopole solutions 
may be represented in the form $ = [F(u + 2r) — F(u)]/r and axisymmetric multipoles may 
be built up by applying the ^-translation operator 

d z = cos 8(d r — d u ) — r _1 sin 8d e (21) 

to these solutions. Then 7 may be obtained via (§). For the calibration measurements in 
Sec. [V], we use the solutions 

$ = (d z y[u(u + 2r)}-\ (22) 

obtained by applying (d z ) e to the fundamental Lorentz invariant solution l/(x a x a ). This 
solution is well behaved above the singular light cone u = 0. 



III. THE NONLINEAR ALGORITHM 

The nonlinear evolution equation (||) can also be recast in the form of (|2"0| ) in terms of 
ip = r^, for an appropriate choice of 2-dimensional wave operator nw. In this case, the 
(u, r) submanifold is not flat so that it would not be appropriate to base □ ( 2 ) upon a flat 
metric. Indeed, such a choice would lead to an improper domain of dependence and could 
not be used as the basis for a stable algorithm. It would seem more natural to use the □ ^ 
operator of the metric induced in the (u,r) submanifold by the 4-dimensional metric (|l]). 
Here we pursue an alternative choice based upon the line element 

da 2 = 2l {lx n v) = e 2f3 du[—du + 2dr], (23) 

where = u tfl is the normal to the outgoing null cones and is the other null vector normal 
to the spheres of constant r. Although this choice is not unique and other possibilities deserve 
exploration, it leads to the simplest 7i-terms when reexpressing @ in the form (|20D . Because 
the domain of dependence of da 2 contains the domain of dependence of the induced metric 
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of the (u, r) submanifold, this approach does not lead to convergence problems associated 
with the CFL condition. 

The wave operator associated with (|23| ) is 



□ {2 V = e~ 2/3 [2^ ru -(^,r), r ] (24) 
and the nonlinear evolution equation (|5|) becomes 

□ cty = e -^n, (25) 

where 

U = "(-).r7 - -P(7,^ + \ sm9(-?-) te )l r 
r r 2 sma 



(j >r U sm9) t9 



_ ri ^p, + lr 3 e 2(^) {Ur) 2 
sin 9 4 

+i e ^-7) [(/3 )2 + sin ( ^) ]. (26 ) 
r sin 9 

Because all 2-dimensional wave operators are conformally flat, with conformal weight —2, 
the surface integral of (^) over a null parallelogram gives an integral equation analogous to 

(0), 

ipQ = + ~ 4>r + 7: / du drH. (27) 

This allows the evolution of 7 by the same basic marching algorithm as described in Sec. |Tl 
for the scalar wave and linearized wave cases. The additional modifications are that /3, U, 
and V must be simultaneously marched out the null cone using the nonlinear hypersurface 
equations (U)-©. Because of the hierarchal structure of these equations, 7, /3, U, and V 
may be marched in that order without any matrix inversions or other implicit operations. 
The basic computational cell and finite difference constructions are described in the next 
section. 

Near the origin, the metric approaches the Minkowski metric so that the stability of the 
nonlinear algorithm in this region should be subject to the same Courant limit as for the 
linearized equations. Near scri, the gravitational variables have the asymptotic form 

7 = K + r- 1 c + 0(r-- 2 ) (28) 
f3 = H + 0{r- 2 ) (29) 
U = L + 0(r- 1 ) (30) 



V 



2 (L sin 1 



sin 9 

+ re W-K> 



sm.9 sm 6 9 



-2e 2H M + 0(r- 1 ), (31) 

where M. corresponds to Bondi's mass aspect. In a standard Bondi frame at scri K, H, 
and L all vanish, but not in null cone coordinates adapted to a Minkowski frame at the 
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origin. This dependence can lead to drastic behavior of the ti-coordinate at large distances. 
In a numerical study of spherically symmetric, self-gravitating scalar radiation fields @, 
H — > oo as a horizon is formed and an infinite redshift develops between central observers 
and observers at scri. In that case, a consideration of the domains of dependence indicates 
that the step size Au for stable evolution must approach zero as the horizon is formed. The 
divergence of the outgoing null cone equals e~ 2f3 jr. If (3 — > oo at a finite value of r then a 
caustic will in general form. When this occurs, a single null cone coordinate system cannot 
cover the entire exterior region of the spacetime. 

In the more general case being considered here, it is also possible for the u-direction to 
become spacelike at large distances, corresponding to the coefficient of du 2 in (|1]) becoming 
negative. (The u = constant hypersurfaces, of course, must remain null.) As discussed in 
Sec. [V], this does not affect the stability of the algorithm. The algorithm is also valid when 
the vertex of the null cone is replaced by an inner boundary consisting of a timelike or null 
worldtube, so that it may also be applied to other versions of the characteristic initial value 
problem ||. 



IV. FINITE DIFFERENCE TECHNIQUES 



The numerical grid is based upon the outgoing null cones, using the compactified radial 
coordinate x = r/(l + r) and the angular coordinate y = — cos9. Thus points at scri are 
included in the grid at x = 1. 

In order to improve numerical accuracy at the grid boundaries, the code is written in 
terms of the variables ip — ry = rj/ sin 2 9, (3, S = (V — r)/r 2 and U = U / sin#. For a pure 
quadrupole mode, 7 has constant angular dependence. 

To develop a discrete evolution algorithm, we work with two sets of spatial grid points, 
both of which have the constant spacing Aa; = 1/N X and Ay = 1/N y . The first grid is 
defined by (u n ,Xi,yj) = (nAu, iAx,jAy). The second grid is shifted (staggered) in both 
the x and y variables and is thus defined by (u n , x i+ ± , = (nAu, (i + \)Ax, (j + |)Ay). 
Note that the staggered grid extends beyond the physical boundary x — 1. This peculiarity 
is successfully exploited for a smooth calculation of the metric at scri. The time step is 
variable and is limited by the largest possible value that would satisfy the CFL condition 
over the entire grid. 

A staggerered grid is not necessary for the scalar wave equation but its introduction 
for the gravitational case is dictated by a detailed von Neumann stability analysis of the 
linearized equations. Accordingly, 7, (3 and S reside on the (xi, yj) grid while U is placed on 
the (x i+ i,yj + i) grid. We denote values of a field F at the site (n, i,j) by F£j = F(u n , Xi, yj). 
We use centered second order differences for derivatives at points not on the edges of the 
grid, e.g. for an arbitrary field F 

F ,y\i,j = 2Ay ~ Fi *-i) ( 32 ) 

To calculate derivatives of the field at the edges of the grid, we use backward and forward 
second-order differences, e.g at the y = ±1 edges of the grid, where j = ±N y 

F ,y\i,±N y = ^^-(-^,±^^2 + 4F it±NyTl - 3F h±Ny ) (33) 
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A. The Hypersurface Equations 

In terms of the numerical variables 7, /3, S and U, and expressed in the coordinates x 
and y used in the code, the hypersurface equations (|2|)-((D are 

(3, x = Hp (34) 
( ^e*V»-«^. 2 _£!_ W! , (35) 

x 2 S,, + -^—S = H s , (36) 
1 — X 

where the source terms 7ig, Tiu and Tis are given by 

H p = l -x{l-x){l-y 2 Y^, x Y (37) 
Hv = P,xy - , , P, y + iy% + (1 - y 2 )[2%((l - y 2 )% - 2y 7 ) - % y ] (38) 

X 1 J. X J 

4 „ * 1 4 - 

ft s = -1 - zy[— -C/ + xU x ] + -(1 - y 2 )x[^— -C/,,, + zE^] 

_i x 4 ( ^ )2(1 _ y2)e 2((l-^)7)-/3 _ e 2( / 3-(l-^ W{ _ 1 _ m _ 2yf3y 

+ (1 - y 2 ) [10 7 + 8^ + 8 7 2 + 4yjP ty + P >yy + 2 

-(1 - y 2 ) 2 (8 7 2 + 2%P tV + j, yy + Syj%) + 2(1 - y 2 f{l,yf}. (39) 

Note that Eq. fl3"4"| ) has just one radial derivative, and can be evaluated at the points (n, i — 
We can discretize it as follows 

A,; = A-ij + (W/3)i~i i (40) 

Equation (^), which contains a second radial derivative, is evaluated at the points 
(n, i — 1, j). Near the origin, where U = 0(x), it becomes a singular differential equation. In 
order to integrate it to second order accuracy it is necessary to apply numerical regularization 
to the derivatives on the left hand side. Noting that d/dx = 4x 3 d/dx 4: , Eq. fl35|) becomes 



2x(l - x)[x 4 U x ]^ + x 2 [l - (1 - x)((3, x - (1 - y 2 )i, x )]U tX 

= e 2 ^ 1 -y 2 *\l-x)H u . (41) 

Using the identity 



x . 



1 - =2Ak ! _ 1 (x _ 1 + 4-1-0, (42) 

2 1 1 2 2 2 



we can discretize (^T|) as 



„2~ T ^2 ^ — \ [4-4 (^*>J _ ~~ 4-1-1 _ ^i-2,i)] 



1 2 1 1 2 

+ -a^Azfl - (1 - ^_i)(Aj - A-2,j - (1 - 2/|)(7i,j - 7i-2j))^](^i,j - ^-2,j 
(Ax) 2 (l-x l _ 1 )(H c/ ),- lj , (43) 



9 



The above is a 3-point formula, and it can not be applied at the points at X{ = Ax, 
however we know the asymptotic behavior of U at the origin, and we can use it to construct 
a starting algorithm for these points. We approximate the Bondi equations for 7 and U by 
the leading two terms in a power series, 

7 = ar 2 + br 3 

U = Ay(ar + \br 2 ), (44) 
Expanding j jU to the same order, we obtain for the rate of change of a and b 

6 h 

a, u = 7° 

5 

b, u = 0. (45) 

By fitting a least square polynomial to 7/r 2 near the origin, we can read off the coefficients 
a and b and evaluate U on the next hypersurface. This approximation is consistent with the 
global second order accuracy of the algorithm. 

As with Eq. (|34]), Eq. (^) can be approximated at sites (n, i — |, j) as follows 

1 ,„ „ , x i-l 



x l 1 "a — — Si^ij) + (Sij + Si-ij) 

l ~2 Ax 1 — x ■ 1 



2 



(«*W (46) 



B. The Evolution Equation 

In practice, the corners of the null parallelogram, P, Q, R and S, cannot be chosen to lie 
exactly on the grid because the velocity of light in terms of the compactified coordinate x is 
not constant even in flat space. Numerical experimentation suggests that a stable algorithm 
with high accuracy results from the choice made in Fig. [l|. The essential feature of this 
placement of the parallelogram with respect to a coordinate cell is that the sides formed 
by incoming rays intersect adjacent u-hypersurfaces at equal but opposite x-displacements 
from the neighboring grid points. Solution for the null geodesies of the metric fl23|) to second 
order accuracy then gives for the coordinates of the vertices 

— Xp = Xr — 

= At*(l - Xi_i) 2 [l + rwCfiZS 1 + S?)/2]/4 

X{ Xq Xg X{ 

= Au(l - Xl ) 2 [l + n{S^ + S? +1 )/2]/4. (47) 

The elementary computational cell consists of the lattice points (n,i,j) and (n, i ± l,j) 
on the "old" hypersurface and the points (n + l,i — 2,j), (n + l,i — l,j) and (n + 
on the "new" hypersurface (and their nearest neighbors in the angular direction). The 
marching algorithm computes the value of the fields at the point (n + in terms of 

their predetermined values at the other points in the cell. 
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The values of ip at the vertices of the parallelogram are approximated to second order 
accuracy by linear interpolation between grid points. Furthermore, cancellations arise be- 
tween these four interpolations so that the evolution equation (j27D is satisfied to fourth order 
accuracy, provided the integral can be calculated to that accuracy. This is accomplished by 
approximating the integrand by its value at the center of the parallelogram. To second order 
accuracy, this gives 

/ dudrH = H c [ dudr = 7- Aw (j-q — r P + r s — r R )H c , (48) 

J A J A ^ 

where the centered value 7i c can be obtained by averaging between appropriate points in 
the cell. Thus the discretized version of (1371) is given by 



C +1 = €-2, C+i, V>I\ C-i) 

+ ^Au (r Q - r P + r s - r R )H c (49) 

where T is a linear function of ip and the j index has been suppressed. 

Consequently, it is possible to move through the interior of the lattice, computing 
explicitly by an orderly radial march. This is achieved by starting at the origin at time 
u n+l . Field values vanish there. Next, proceed outward one radial step using the boundary 
conditions (discussed below). Then step outward to the next interior radial point using (fl9"|), 
iterating this process throughout the interior and for all angles. This updates all field values 
stretching to scri along the new null cone at u n+1 , thus completing one evolutionary time 
step. 

The above scheme is sufficient for accurate evolution in a neighborhood of the origin. 



Global evolution, including the points at x — 1, requires careful manipulation of (^) to 
avoid problems from the fact that ip ~ rK at scri. Thus the direct use of this formula is 
not possible for the point at scri, while points near scri would suffer serious loss of accuracy. 



We renormalize ([4"9"D in the following way. First, we introduce the quantity (f) = — x) 



Near scri has the desired finite behavior, while near the origin it leaves unchanged the 
constant coefficient form of the evolution equation, thus preserving the stability properties. 
With this substitution and with the use of (f48"D, the evolution equation (^) becomes 



(f) Q = \auxqH c + j- — ^j-(0p - \&ux P TL c ) 

4 (1 — Xp) 4 

+ Q ~ {<Ps + \auxs7Q ~ n I l 9 l ftto + \^ux R H c ) (50) 

Now all terms have finite asymptotic value. The coefficient (1—xq) / (1—xs) has 0/0 behavior 
at scri but approaches the limit 1. Further refinement is possible with the use of the explicit 
second order approximation for the characteristics ( |4~7|) which leads to the approximation 

( / ~ ^ = 1 + 2— g— (51) 
(l-x s ) 1-6 1 ; 

where 
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6= 1 -Au(l + x t (S?+ 1 1 + S? +1 )/2) (52) 

The final result is that the equation (|50"D propagates <fi radially outward one cell with an 
error of fourth order in grid size. This is valid for all interior points and the point at scri. 
The error in each cell compounds to a third order error on each null cone and a second order 
global error after evolving for a given physical time. Second order global accuracy is indeed 
confirmed by the convergence tests described in Sec. |V|. 

We mentioned that a modified form of the basic grid cell Eq. fl4~7|) is used at the origin. 
This is necessary since the incoming characteristic through the points P and R can not be 
centered at x — 0. The corners of the modified cell are given by 

x P = 
1 . 

Xr = - Au 

xi - xq = x s - xi = ^A«(l - Xi) 2 . (53) 

Only the linear terms of 7i are kept while evolving the first point, i.e. for x\ = Ax. This 
reduces equation (|27D to 

If 1 

i>Q = 4>s ~ 4>r ~ j J dudr-(r 2 U, y ), r . (54) 



Using the expansion ([P]) for U near the origin, the integral simplifies further to 

f 12 

/ dudr(3ar-\ br 2 ). (55) 

J a 5 

The integrand is now evaluated to second order accuracy at u n+ ^ = u n + 4p using 



.1 6 Au 

2 = a + -b 

5 2 



i 



and b n+ 2 = b n . Keeping higher order terms would not improve the global convergence rate 
of the code. 



V. CODE TESTS 
A. Testbeds 

As we have shown in Sec. U linearized solutions of the Bondi equations can be generated 
by solutions of the scalar wave equation, thus supplying a complete set of test beds for the 



very weak regime. For the nonlinear case, exact boost and rotation symmetric solutions [10 



of the Bondi initial hypersurface equations have also been found UlTJl . They have been used 
to check the radial integrations leading from 7 to (3, U and V but they do not provide a 
test of the evolution algorithm. However, in the course of this work, we have found that 
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one of these initial data sets is in fact preserved under time evolution and is an exact static 
solution of the nonlinear vacuum equations. 

This solution provides an important test bed for null cone evolution codes. Except for 
spherically symmetric cases, it is the only known solution of Einstein's equation which can 
be expressed explicitly in null cone coordinates with no singularity at the vertex. It has the 
form 

2 e 7 = 1 + E 

(1 + S) 2 



e 2 ? 



U 



4S 



a 2 p \fr 2 — p 2 



r S 

y = r (2a 2 p 2 -a 2 r 2 + l) ^ 



where p = rsinO, S = y 1 + a 2 p 2 and a is a free scale parameter. It is remarkable that 
the null data 7 is time independent under evolution, as can be verified by inserting the 
above expressions in (||). Because this solution is static, as well as boost and rotation 
symmetric, the commutator between the boost and time translation symmetries implies 
that it has an additional translation symmetry in the boost direction. Thus the solution 
falls into several overlapping and widely studied classes, including the static cylindrically 
symmetric spacetimes and the Weyl static axisymmetric spacetimes. However, this solution, 
which we call SIMPLE, has not previously been singled out, apparently because it cannot 
easily be identified in the traditional coordinates used for studying static solutions. Because 
of its cylindrical symmetry, it is clear that this solution is not asymptotically flat but it 
can be used to construct an asymptotically flat, nonsingular solution by smoothly pasting 
asymptotically flat null data to it outside some radius R. The resulting solution will be 
static and given by (|56|) in the domain of dependence interior to R. Numerical solutions 
generated by this technique are used in the code calibration tests presented below. 

In addition, global energy conservation provides an important test bed. The Bondi mass 
loss formula is not one of the equations used in the evolution algorithm but follows from 
those equations as a consequence of a global integration of the Bianchi identities. 



B. Convergence and Stability 

We have tested the algorithm to be second order accurate and stable, subject to the CFL 
condition, throughout the regime in which caustics and horizons do not form. In Sec. |T|, we 
showed how the the linearized Bondi equations may be reduced to the scalar wave equation 
by local operations. For very weak data, the nonlinear equations approximate the linear 
equations so that we would expect the global stability of the nonlinear algorithm to be 
related to the CFL condition for the scalar wave algorithm. Near the origin, stability checks 
show that the time step is limited by (|19|) with K — 8, which is twice the limit found for 
the scalar wave algorithm. This factor of two apparently arises from the use of a staggered 
grid in the gravitational case, which effectively doubles the value of r at which the main 
algorithm takes over from the start up algorithm at the origin. This gives some reassurance 
that the scalar wave algorithm has been optimally adapted to the Bondi equations. 
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By construction, the w-direction is timelike at the origin where it coincides with the 
worldline traced out by the vertex of the outgoing null cone. But even for weak fields, the it- 
direction becomes spacelike at large distances along a typical outgoing ray. This can be seen 
from the metric coefficient g uu = (V/r)e 2f3 — U 2 r 2 e 21 which at large r becomes dominated by 
the the asymptotic behavior U = L + 0(l/r). Geometrically, this reflects the property that 
scri is itself a null hypersurface so that all internal directions are spacelike, except for the null 
generator. For a flat space time, the u-direction picked out at the origin corresponds to the 
null direction at scri but it becomes asymptotically spacelike under the slightest deviation 
from spherical symmetry. 

By choosing initial data of very small amplitude (|7| ~ 10 -9 ), we have performed conver- 
gence tests of the numerical solutions against the solutions of the linearized equations. The 
linearized solutions (^) were given as initial data at u = and we compared the numerically 
evolved solutions to the linearized solutions at a central time of u = 0.5. We observed that 
for the low angular momentum solutions (£ =2, 3, 4) the code is superaccurate, i.e. the solu- 
tions converge to the exact result at a rate faster than second order in the grid size. This is 
to be expected since for these solutions the hatted variables used in the code exhibit angular 
dependence that is at most quadratic in y, so that the second order accurate y-derivatives 
are calculated exactly. The error, as measured by the L 2 norm, of a numerical solution with 
higher harmonics (£ = 6) is graphed in Fig. |2|. The slope of the graph gives a convergence 
rate of 2.04 ± 0.01 with respect to grid size. This result is insensitive to the particular norm 
used, i.e. we also verified second order convergence in the L\ and norms. 

Second order convergence has also been checked against the exact static solution SIM- 
PLE. Since this solution is not asymptotically flat, we match the initial data smoothly to 
asymptotically flat data for a nonstatic exterior. Consideration of the domain of dependence 
implies that the matching boundary propagates along an ingoing null hypersurface. Thus, 
we can obtain a reliable measure of how accurately the evolution preserves the static inte- 
rior, provided we restrict the calculation of the error norm to a region not yet influenced by 
the exterior nonstatic data. We matched one such static solution in the interior (x < 0.5) 
to smooth exterior data with compact support. We calculated the L ro norm for the region 
x < 0.4 at time u=0.25, and considered the dependence of the error on the grid size. The 
preservation of the static interior to graphical accuracy is shown in Fig. |3|, while the second 
order convergence of the error is demonstrated in Fig. |4j. In addition, for a wide variety of 
initial data having unknown analytic solution, we have verified that the numerical solution 
converges to second order in the sense of Cauchy convergence. 

Stability of the code in the low to medium amplitude regime has been verified experimen- 
tally by running arbitrary initial data until it radiates away to scri. At higher amplitudes, 
it is expected that physical singularities will arise, but we have not yet explored this regime. 
Figure [5] shows a sequence of time slices of the numerical evolution of some arbitrary initial 
data of compact support. Note the rich angular structure that arises at u « 0.25 and then 
dissipates. At u ~ 1.5 the amplitude of the field is sufficiently small so that it appears to 
be zero in the figure. It continues to decay at later times. 
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C. Energy Conservation 



The Bondi mass loss formula relates the gravitational radiation power to the square of 
the news function. It follows from the equations used in the algorithm as a consequence 
of a global integration of the Bianchi identities. Thus it not only furnishes a valuable tool 
for physical interpretation but it also provides a very important calibration of numerical 
accuracy and consistency. 

Historically, numerical calculations of the Bondi mass M B have been frustrated by tech- 
nical difficulties arising from the necessity to pick off nonleading terms in an asymptotic 
expansion about infinity. For example, the mass aspect A4 must be picked off in the the 
asymptotic expansion ( [31]) for V. This is similar to the experimental task of determining the 
mass of an object by measuring its far field. In the non-radiative case it can be accomplished 
by measuring gravity gradients, but otherwise this approach can be swamped by radiation 
fields. In the computational problem, further complications arise from gauge terms which 
dominate asymptotically even over the radiation terms. We have recently developed a second 



order accurate algorithm for calculating the Bondi mass [ 12j . It avoids the above problems 
through the use of Penrose compactification and the introduction of renormalized variables 
in which Bondi's mass aspect appears as the leading asymptotic term. The Bondi mass 
algorithm depends only upon fields on a single null hypersurface. It has been incorporated 
into the present evolution code to calculate the mass at any given retarded time. 
In the present formalism, the news function iV is given by [|TJ 



2e 2H N = 2c u + 



(sinflc 2 L), 
sin 9 c 



+ e~ 2K uj sin 9 



sin 9 ijj 2 



(57) 



Here u is the conformal factor relating the asymptotic 2-geometry to the unit sphere geom- 
etry of a Bondi frame, i.e. 



e 2K d9 2 + sin 2 



-2K ji2 



uo- 2 {d9 2 B + sm 2 9 B d<p 2 B ), 



(58) 



where 9b and <ps = 4> are Bondi spherical coordinates. Calculation of u complicates the 
calculation of the news function. The simplest approach is to set y = — cos 9 and yg = 
— cos 9b- Then 
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dy_ 



B 



where 



ys = tanh 



This gives 



dy 

y dy 



,2if 



o 



(59) 



(60) 



2e 



K 



[l + y)e* + (l-y)e-* 



(61) 



where 
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A-/*^- (62) 

In order to prepare this integral in an explicitly regular form for computation we introduce 
an auxiliary parameter a and rewrite (|62l) as the double integral 



A = 2 f dy I da e 2aK K, (63) 
Jo Jo 

where K = Kj (1 — y 2 ) is regular at the poles. It is then straightforward to obtain a second 
order accurate finite difference formula for the news function. 

The Bondi formula for energy conservation between central times Mo and u takes the 
form C =0, where 



C = M B (u) - M B (u ) + \ f dy f due^u^N 2 . 

2 J-l Ju 



(64) 



Figure [| graphs C relative to the initial mass Mb(uq) for a numerical evolution of the 
polynomial data 



[(x - x l )(x-x 2 )(y 2 - ?/o)] f 
[(xi - x 2 )yo\ 



7 = A ~ r^T2 ( 65 ) 



with compact support in the domain (xi < x < x 2 ) X (— yo < y < yo) and amplitude 
parameter A. For the graph, we have chosen A = 0.3, x\ = 0.1, x 2 = 0.5 and y = 0.5, and 
evolved the numerical solution up to u = 0.01, on a grid of 512 radial x 128 angular points. 
The bulk of the error occurs in the calculation of the Bondi mass, whose accuracy is more 
sensitive to grid size than the accuracy of either the news function or the evolution code. 
Most of this error may be removed by using Richardson extrapolation to take advantage 
of the known second order accuracy of the Bondi mass. For example, if F n (xi) is a second 
order accurate finite difference approximation to the function f(x) on a grid of n points, 
then (4F n — F n / 2 )/3 approximates / to third order and in fact to fourth order if odd orders 
are absent in the approximation. This absence of odd orders indeed holds for the Bondi 
mass because all derivatives, interpolations and integrals are centered. Thus, introduction 
of subgrids obtained by subsampling, leads to a fourth order expression for M B , with the 
corresponding relative error in energy conservation also graphed in Fig. ^. In this way, energy 
conservation is attained to 0.4% accuracy. Note that only a single evolution on a fixed grid 
is necessary here because Richardson extrapolation is applied when calculating the Bondi 
mass. For the purposes of the graph, we have done this for each time that C is plotted, but 
to check energy conservation, it suffices to do it only at the initial and final times. Figure |6| 
serves as a rewarding testament to the virtues of a code with known convergence rates. 



VI. CONCLUSION 

We have constructed a second order accurate evolution algorithm for the null cone initial 
value problem for axisymmetric vacuum spacetimes. Energy conservation is maintained to 
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second order accuracy. Extensive tests of the algorithm establish that it is globally valid in 
the regime where horizons and caustics do not develop. This generates a large complement of 
highly accurate numerical solutions for the class of asymptotically flat, axisymmetric vacuum 
spacetimes, for which no analytic solutions are known. All results of numerical evolutions 
in this regime are consistent with the theorem of Christodoulou and Klainerman that 
weak initial data evolve asymptotically to Minkowski space at late time. The code is now 
being tested in the strong field regime for application to the study of black hole formation. 
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APPENDIX: 

We sketch here the von Neumann stability analysis of the algorithm for the linearized 
Bondi equations. The analysis is based up freezing the explicit functions of r and y that 
appear in the equations, so that it is only valid locally for grid sizes satisfying Ar « r and 
Ay << 1. However, as is usually the case, the results are quite indicative of the stability of 
the actual global behavior of the code. 

Starting with the hatted code variables introduced in Sec. [IV] and setting T = r 2 U and 



G = ry, the linearized Bondi equations ([]) and (|7|) take the form 

r 2 r rr -2T = 2[4y - (1 - y 2 )d y ](rG, r - G) (Al) 

and 

2G >ur - G, rr = -(l/2r)r ry . (A2) 

Freezing the explicit factors of r and y at r = R and y = Y, introducing the Fourier modes 
G = e su e lkr e l[ y (with real k and I) and setting V = AG, these equations imply 

A = 2(1 - ikR)[4Y - (1 - Y 2 )il]/{2 + R 2 k 2 ) (A3) 

and 

Ais = -2k + Al/R. (A4) 

For stable modes, Re(s) > 0. This requires that the llm(A) < which will be satisfied unless 
Ykl < 0. In the latter case, unstable solutions exist to the PDE's obtained by freezing the 
coefficients in the linearized equations (|A1| ) and ( |A2| ). The linearized equations themselves 



do not have unstable modes but they arise in the frozen coefficient formalism from dropping 
the boundary condition of spherical topology on the y-dependence. For a global solution, 
G should not have periodic dependence on y but instead be decomposed into spin-weight 
2 harmonics, in which case instabilities would not arise in the above analysis. Thus these 
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unstable modes of the frozen PDE are artificial and should be discarded by requiring Ykl > 
when analyzing the stability of the corresponding FDE. 

Consider now the FDE obtained by putting G on the grid points rj and T on the staggered 
points r I+ i/ 2 , while using the same angular grid yj and time grid U]y. Let P, Q, R and S 
be the corner points of the null parallelogram algorithm, placed so that P and Q are at 
level N + 1, R and S are at level N, and so that the line PR is centered about rj and 
QS is centered about r/ + i. For simplicity, we display the analysis at the equator where 
Y — 0. Then, using linear interpolation and centered derivatives and integrals, the null 
parallelogram algorithm for the frozen version of the linearized equations leads to the FDE's 

(R/Ar) 2 (T I+3/2 - 2F I+1/2 + r 7 _ 1/2 ) - (F I+3/2 + V^ l/2 ) 

= -6 y [2(R/Ar)(G I+1 - Gj) - (G J+1 + G T )\ (A5) 

(all at the same time level) and 

W+l ~~ W ~~ W+l + °7 

+ (An/4Ar)(-Gf+ 1 + 2G? +1 - Gf + 1 - Gf +2 + 2G? +1 - G?) 

= -(Au/sRfa&Wp - rf+) 2 + rf +3/2 - rf +1/2 ), (A6) 

where S y represents a centered first derivative. Again setting T = AG and introducing the 
discretized Fourier mode G = e sN ^ u e lkI ^ r e lU ^v ^ we have 5 y = i sin (I Ay) /Ay and (|A5|) and 
( |A6| ) reduce to 

A[(R/Arf(l - cos a) + cos a] = -L[{2R/Ar) sin(a/2) + i cos(a/2)] (A7) 

and 

e sAu = _ e ia( C * _ up _ ( A g) 

where L = sm(lAy) j Ay, a = kAr, C = ie ia l 2 sin(a/2) + (Aw/4Ar)(l — cosa) and D = 
(LAu/8R) sin(a/2). The stability condition that Re(s) < then reduces to Re[CD(A — 
A*)} > which is equivalent to 1 + cosafl — (Ar/R) 2 } > 0. Thus this stability condition is 
automatically satisfied and poses no constraint on the algorithm. 

The corresponding analysis at the poles Y = ±1 again leads to (JA8|), where now 

A[(R/Ar) 2 (l - cos a) + cos a] = -4Y[2i(R/Ar) sin(a/2) - cos(a/2)]. (A9) 

The stability condition Re[CD(A — A*)] > is satisfied provided Ykl > 0, which rules out 
the artificially unstable solutions of the frozen PDE discussed above. 

As a result, local stability analysis places no constraints on the algorithm. This may seem 
surprising because not even the analogue of a CFL condition on the time step arises but it 
can be understood in the following vein. The local structure of the code is implicit, since it 
involves 3 points at the upper time level. Implicit algorithms do not necessarily lead to a CFL 
condition. However, the algorithm is globally explicit in the way that evolution proceeds by 
an outward radial march from the origin. It is this feature that necessitates a CFL condition 
in order to make the numerical and physical domains of dependence consistent. 
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FIGURES 




FIG. 1. Line segments drawn at forty- five degrees represent radial characteristics. Their in- 
tersection defines the fundamental null parallelogram PQRS shown superimposed upon the com- 
putational cell, which consists of the points marked by circles and their nearest neighbors in the 
angular direction. 
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FIG. 2. The error, as measured by the L2 norm, of a numerical solution with higher harmonics 
(£ = 6). The computation is made on grids of size N x equal to 24, 48, 72, 96 and 120, while keeping 
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FIG. 3. Evolution of initial data given by SIMPLE in the interior region, and patched smoothly 
to an asymptotically flat exterior. The static interior is preserved to graphical accuracy. 



22 



-5.5 




-7.5 1 1 1 1 1 1 1 1 1 1 1 

-2.6 -2.4 -2.2 -2.0 -1.8 -1.6 

log 10 (Ax) 



FIG. 4. The error in the evolution of the initial data of Fig. |3| up to u = 0.25, as measured 
by the norm. The error is computed on grids of size N y equal to 16, 24, 32, 48 and 64, 
while keeping N x = 3N y . The convergence rate is 1.92, in good agreement with the theoretically 
expectation of second order accuracy. 



23 



11=0.0000 



11=0.2665 




(d) 

u=1.0660 u=1.4657 




FIG. 5. A sequence of time slices of the numerical evolution of initial data of compact support. 
Note all the angular structure that arises at about u = 0.25, which later decays. At u « 1.5 the 
amplitude of the field has decayed below those values which can be observed in the figure. 
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FIG. 6. Graph of the relative error in C calculated up to u = 0.01 for a numerical evolution of 
the data (|65|), on a grid of 512 radial x 128 angular points. The circles show the error as calculated 
from the computed values of the mass, while the squares show the error after using Richardson 
extrapolation, based on the known convergence rate of the algorithm. 
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